SPARS trial B
Descriptive plots of the stimulus-response relationship for the SPARS, NRS_NP, and CNRS_P
Peter Kamerman and Tory Madden
14 February 2018
This script is part 1 of our analysis of the stimulus-response characteristics of the SPARS, NRS_NP, and CNRS_P. This script generates exploratory plots of the relationship between stimulus intensity and scale ratings.
The three scales measure the following ranges of somatic sensation:
CNRS_P: 0 (no pain) to 100 (worst pain you can imagine)
NRS_NP: 0 (no sensation) to 100 (pain)
SPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)
Because the the stimulus range was calibrated to the sensitivity of each participant (compared to the fixed range of intensities used in Trial A), all analyses use the rank order of the nine stimulus intensities each participant was exposed to rather than the absolute intensities of the stimuli used.
For convenience, experimental blocks have been labelled in the order they took place within each scale. Each scale was tested across four successive blocks of 27 stimuli each. The order of in which the scales were assessed was randomized.
Import and inspect data
# Import
data <- read_rds('./data-cleaned/SPARS_B.rds')
# Inspect
glimpse(data)## Observations: 2,256
## Variables: 9
## $ PID <chr> "ID01", "ID01", "ID01", "ID01", "ID01", "ID01"...
## $ scale <chr> "SPARS", "SPARS", "SPARS", "SPARS", "SPARS", "...
## $ block_number <int> 2, 2, 2, 4, 4, 4, 6, 6, 6, 8, 8, 8, 11, 11, 11...
## $ trial_number <int> 9, 15, 23, 7, 20, 25, 9, 18, 22, 3, 17, 23, 3,...
## $ intensity <dbl> 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25...
## $ intensity_char <chr> "2.25", "2.25", "2.25", "2.25", "2.25", "2.25"...
## $ intensity_rank <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ rating <int> -31, -20, -48, -48, -21, -23, -48, -45, -47, -...
## $ rating_positive <dbl> 19, 30, 2, 2, 29, 27, 2, 5, 3, 0, 1, 3, 50, 50...Clean and transform data
############################################################
# #
# Clean #
# #
############################################################
data %<>%
# Select required columns
select(PID, scale, block_number, intensity, intensity_char,
intensity_rank, rating, rating_positive)
############################################################
# #
# Calculate 'Tukey trimean' #
# #
############################################################
# Define tri.mean function
tri.mean <- function(x) {
# Calculate quantiles
q1 <- quantile(x, probs = 0.25, na.rm = TRUE)[[1]]
q2 <- median(x, na.rm = TRUE)
q3 <- quantile(x, probs = 0.75, na.rm = TRUE)[[1]]
# Calculate trimean
tm <- (q2 + ((q1 + q3) / 2)) / 2
# Convert to integer
tm <- as.integer(round(tm))
return(tm)
}
############################################################
# #
# Generate core data #
# #
############################################################
# Add sequential block numbers, tagged by scale
data %<>%
# Group by ID and scale
group_by(PID, scale) %>%
# Arrange blocks
arrange(block_number) %>%
# Add sequential
mutate(block_sequential = dense_rank(block_number),
block_sequential = paste0(scale, ' - ', block_sequential)) %>%
ungroup()
# Convert NRS_NP and CNRS_P to SPARS-equivalent ranges (+50 to 50)
## Please see section: Equivalent unit rating for details
data %<>%
# Fix column format to do the maths
mutate(rating = as.numeric(rating)) %>%
# CNRS_P conversion(divide rating by 2)
# NRS_NP conversion (-100 and divide then by 2)
mutate(rating_equivalent = case_when(
scale == 'CNRS_P' ~ rating / 2,
scale == 'NRS_NP' ~ (rating - 100) / 2,
TRUE ~ rating
))
# Calculate the participant trimean for each scale (using 'rating_equivalent')
data_tm <- data %>%
group_by(PID, scale, intensity_rank) %>%
summarise(tri_mean = tri.mean(rating_equivalent)) %>%
ungroup()Exploratory plots
Raw intensity ratings
Data are shown on their original scales: CNRS_P: 0 (no pain) to 50 (worst pain you can imagine), NRS_NP: -50 (no sensation) to 0 (pain), and SPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine).
# Generate plots
plot_raw <- data %>%
group_by(PID) %>%
nest() %>%
mutate(plots = map2(.x = data,
.y = PID,
~ ggplot(data = .x) +
aes(x = intensity_rank,
y = rating,
colour = scale,
fill = scale) +
geom_smooth(se = FALSE) +
geom_point() +
scale_y_continuous(limits = c(-50, 100)) +
scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
scale_fill_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
scale_colour_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
labs(title = paste0(.y, ': Raw participant-level stimulus-response ratings on the CNRS_P, NRS_NP and SPARS'),
subtitle = 'Plots are conditioned on experimental block\nCNRS_P: 0 (no pain) to 100 (worst pain you can imagine)\nNRS_NP: 0 (no sensation) to 100 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
x = 'Stimulus intensity (rank)',
y = 'Intensity rating') +
geom_hline(yintercept = 0, linetype = 2) +
facet_wrap(~ block_sequential)))
# Print plots
walk(plot_raw$plots, ~ print(.x))
Equivalent units ratings
Raw scores for the CNRS_P and NRS_NP scales (both rated on 0 to 100 scales) were converted to SPARS equivalent ranges (-50 to +50). The scaling of the CNRS_P and NRS_NP were as follows:
CNRS_P: raw 0 to 100 scores were converted to a 0 to 50 range by dividing 2, such that after the scaling, 0 = no pain and 50 = worst pain you can imagine. This is the positive range of the SPARS.
NRS_NP: raw 0 to 100 scores were converted to a -50 to 0 range by subtracting 100, and then dividing 2, such that after the scaling, -50 = no sensation and 0 = pain. This is the negative range of the SPARS.
# Generate plots
plot_equi <- data %>%
group_by(PID) %>%
nest() %>%
mutate(plots = map2(.x = data,
.y = PID,
~ ggplot(data = .x) +
aes(x = intensity_rank,
y = rating_equivalent,
colour = scale,
fill = scale) +
geom_smooth(se = FALSE) +
geom_point() +
scale_y_continuous(limits = c(-50, 50)) +
scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
scale_fill_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
scale_colour_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
labs(title = paste0(.y, ': Scaled participant-level stimulus-response ratings on the CNRS_P, NRS_NP and SPARS'),
subtitle = 'Plots are conditioned on experimental block\nCNRS_P: 0 (no pain) to 50 (worst pain you can imagine)\nNRS_NP: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
x = 'Stimulus intensity (rank)',
y = 'Scaled intensity rating (-50 to 50)') +
geom_hline(yintercept = 0, linetype = 2) +
facet_wrap(~ block_sequential)))
# Print plots
walk(plot_equi$plots, ~ print(.x))
Tukey trimean plots
# Generate plot
plot_tm <- data_tm %>%
# Nest
group_by(PID) %>%
nest() %>%
# Plot
mutate(plot = map2(.x = data,
.y = PID,
~ ggplot(data = ..1) +
aes(x = intensity_rank,
y = tri_mean) +
geom_smooth(method = 'loess',
se = FALSE,
colour = '#666666') +
geom_point(shape = 21,
size = 3,
fill = 'orange') +
scale_fill_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
scale_colour_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
scale_y_continuous(limits = c(-50, 50)) +
labs(title = paste0(.y, ': Scaled participant-level stimulus-response rating Tukey trimeans on the CNRS_P, NRS_NP and SPARS'),
subtitle = 'Plots are conditioned on the three scales\nCNRS_P: 0 (no pain) to 50 (worst pain you can imagine)\nNRS_NP: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
x = 'Stimulus intensity (rank)',
y = 'Scaled intensity rating (-50 to 50)') +
geom_hline(yintercept = 0, linetype = 2) +
facet_wrap(~ scale, ncol = 3)))
plot_tm$plot[[2]]
# Print plots
walk(plot_tm$plot, ~ print(.x))
Summary plot (equivalent units)
# Generate plot
data %>%
ggplot(data = .) +
aes(x = factor(intensity_rank),
y = rating_equivalent,
colour = scale,
fill = scale) +
geom_boxplot(alpha = 0.6) +
scale_fill_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
scale_colour_brewer(name = 'Scale',
type = 'qual',
palette = 'Dark2') +
labs(title = 'Scaled group-level stimulus-response ratings on the CNRS_P, NRS_NP and SPARS',
subtitle = 'Plots are conditioned on experimental block\nCNRS_P: 0 (no pain) to 50 (worst pain you can imagine)\nNRS_NP: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
x = 'Stimulus intensity (rank)',
y = 'Scaled intensity rating (-50 to 50)') +
geom_hline(yintercept = 0, linetype = 2)
Session information
sessionInfo()## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2 forcats_0.2.0 stringr_1.2.0
## [4] dplyr_0.7.4 purrr_0.2.4 readr_1.1.1
## [7] tidyr_0.8.0 tibble_1.4.2 ggplot2_2.2.1.9000
## [10] tidyverse_1.2.1 magrittr_1.5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.15 RColorBrewer_1.1-2 cellranger_1.1.0
## [4] pillar_1.1.0 compiler_3.4.3 plyr_1.8.4
## [7] bindr_0.1 tools_3.4.3 digest_0.6.15
## [10] lubridate_1.7.1 jsonlite_1.5 evaluate_0.10.1
## [13] nlme_3.1-131 gtable_0.2.0 lattice_0.20-35
## [16] pkgconfig_2.0.1 rlang_0.1.6 psych_1.7.8
## [19] cli_1.0.0 rstudioapi_0.7 yaml_2.1.16
## [22] parallel_3.4.3 haven_1.1.1 xml2_1.2.0
## [25] httr_1.3.1 knitr_1.19 hms_0.4.1
## [28] tidyselect_0.2.3 rprojroot_1.3-2 grid_3.4.3
## [31] glue_1.2.0 R6_2.2.2 readxl_1.0.0
## [34] foreign_0.8-69 rmarkdown_1.8 modelr_0.1.1
## [37] reshape2_1.4.3 backports_1.1.2 scales_0.5.0.9000
## [40] htmltools_0.3.6 rvest_0.3.2 assertthat_0.2.0
## [43] mnormt_1.5-5 colorspace_1.3-2 labeling_0.3
## [46] stringi_1.1.6 lazyeval_0.2.1 munsell_0.4.3
## [49] broom_0.4.3 crayon_1.3.4